Открытый курс по машинному обучению. Сессия № 2

</center> Автор материала: программист-исследователь Mail.ru Group, старший преподаватель Факультета Компьютерных Наук ВШЭ Юрий Кашницкий. Материал распространяется на условиях лицензии Creative Commons CC BY-NC-SA 4.0. Можно использовать в любых целях (редактировать, поправлять и брать за основу), кроме коммерческих, но с обязательным упоминанием автора материала.

Домашнее задание № 3

Деревья решений для классификации и регрессии

В этом задании мы разберемся с тем, как работает дерево решений в задаче регрессии, а также построим (и настроим) классифицирующие деревья решений в задаче прогнозирования сердечно-сосудистых заболеваний. Заполните код в клетках (где написано "Ваш код здесь") и ответьте на вопросы в веб-форме.


In [2]:
import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
%matplotlib inline
from sklearn.model_selection import train_test_split, GridSearchCV, cross_val_score
from sklearn.metrics import accuracy_score
from sklearn.tree import DecisionTreeClassifier, export_graphviz

1. Простой пример восстановления регрессии с помощью дерева решений

Рассмотрим следующую одномерную задачу восстановления регрессии. Неформально, надо построить функцию $a(x)$, приближающую искомую зависимость $y = f(x)$ в терминах среднеквадратичной ошибки: $min \sum_i {(a(x_i) - f(x_i))}^2$. Подробно мы рассмотрим эту задачу в следующий раз (4-я статья курса), а пока поговорим о том, как решать эту задачу с помощью дерева решений. Предварительно прочитайте небольшой раздел "Дерево решений в задаче регрессии" 3-ей статьи курса.


In [2]:
X = np.linspace(-2, 2, 7)
y = X ** 3

plt.scatter(X, y)
plt.xlabel(r'$x$')
plt.ylabel(r'$y$');


Проделаем несколько шагов в построении дерева решений. Исходя из соображений симметрии, выберем пороги для разбиения равными соответственно 0, 1.5 и -1.5. Напомним, что в случае задачи восстановления регрессии листовая вершина выдает среднее значение ответа по всем объектам обучающей выборки, попавшим в эту вершину.

Итак, начнём. Дерево глубины 0 состоит из одного корня, который содержит всю обучающую выборку. Как будут выглядеть предсказания данного дерева для $x \in [-2, 2]$? Постройте соответствующий график.


In [3]:
plt.scatter(X, y)

X_l = X[X < 0]
y_l = [y[X < 0].mean()] * X_l.shape[0]

plt.plot(X, [y.mean()] * X.shape[0])

plt.xlabel(r'$x$')
plt.ylabel(r'$y$')


Out[3]:
<matplotlib.text.Text at 0x7f8d3ddd96d8>

Произведем первое разбиение выборки по предикату $[x < 0]$. Получим дерево глубины 1 с двумя листьями. Постройте аналогичный график предсказаний для этого дерева.


In [4]:
plt.scatter(X, y)

X_l = X[X < 0]
y_l = [y[X < 0].mean()] * X_l.shape[0]


X_r = X[X >= 0]
y_r = [y[X >= 0].mean()] * X_r.shape[0]

plt.plot(X_l, y_l)
plt.plot(X_r, y_r)

plt.xlabel(r'$x$')
plt.ylabel(r'$y$')


Out[4]:
<matplotlib.text.Text at 0x7f8d3ddd9128>

В алгоритме построения дерева решений признак и значение порога, по которым происходит разбиение выборки, выбираются исходя из некоторого критерия. Для регрессии обычно используется дисперсионный критерий: $$Q(X, j, t) = D(X) - \dfrac{|X_l|}{|X|} D(X_l) - \dfrac{|X_r|}{|X|} D(X_r),$$ где $X$ – выборка, находящаяся в текущей вершине, $X_l$ и $X_r$ – разбиение выборки $X$ на две части по предикату $[x_j < t]$ (то есть по $j$-ому признаку и порогу $t$), а $D(X)$ – дисперсия ответов на выборке $X$: $$D(X) = \dfrac{1}{|X|} \sum_{x_j \in X}(y_j – \dfrac{1}{|X|}\sum_{x_i \in X}y_i)^2,$$ где $y_i = y(x_i)$ – ответ на объекте $x_i$. При каждом разбиении вершины выбираются признак $j$ и значение порога $t$, максимизирующие значение функционала $Q(X, j, t)$.

В нашем случае признак всего один, поэтому $Q$ зависит только от значения порога $t$ (и ответов выборки в данной вершине).

Постройте график функции $Q(X, t)$ в корне в зависимости от значения порога $t$ на отрезке $[-1.9, 1.9]$.


In [5]:
def dispersion(X, y):
    r = 0
    i = 0
    c = 1 / X.shape[0]         
    return c * sum([(_y - c * sum(y)) ** 2 for _y in y])
    

def regression_var_criterion(X, y, t):
    X_l = X[X < t]
    y_l = y[X < t]
    X_r = X[X >= t]
    y_r = y[X >= t]
    
    d = dispersion(X, y) - (X_l.shape[0] / X.shape[0] * dispersion(X_l, y_l)) - (X_r.shape[0] / X.shape[0] * dispersion(X_r, y_r))
    
    return d

In [6]:
t = np.arange(-1.9, 1.9, 0.1)
criteria = [regression_var_criterion(X, y, _t) for _t in t]

plt.plot(t, criteria)

plt.xlabel(r'$x$')
plt.ylabel(r'$y$')


Out[6]:
<matplotlib.text.Text at 0x7f8d3dd6b6a0>

Вопрос 1. Оптимально ли с точки зрения дисперсионного критерия выбранное нами значение порога $t = 0$?

  • Да
  • Нет

Теперь произведем разбиение в каждой из листовых вершин. В левой (соответствующей ветви $x < 0$) – по предикату $[x < -1.5]$, а в правой (соответствующей ветви $x \geqslant 0$) – по предикату $[x < 1.5]$. Получится дерево глубины 2 с 7 вершинами и 4 листьями. Постройте график предсказаний этого дерева для $x \in [-2, 2]$.


In [7]:
plt.scatter(X, y)

print(type(X))

X_l = X[X < 0]
y_l = [y[X < 0].mean()] * X_l.shape[0]

X_ll = X[X < -1.5]
y_ll = [y[X < -1.5].mean()] * X_ll.shape[0]

X_lr = X[(X >= -1.5) & (X < 0)]
y_lr = [y[(X >= -1.5) & (X < 0)].mean()] * X_lr.shape[0]

X_r = X[X >= 0]
y_r = [y[X >= 0].mean()] * X_r.shape[0]

X_rl = X[(X >= 0) & (X < 1.5)]
y_rl = [y[(X >= 0) & (X 
                      < 1.5)].mean()] * X_rl.shape[0]

X_rr = X[(X >= 1.5)]
y_rr = [y[(X >= 1.5)].mean()] * X_rr.shape[0]

X_ = np.r_[X_ll, X_lr, X_rl, X_rr]
y_ = np.r_[y_ll, y_lr, y_rl, y_rr]

plt.plot(X_, y_)

plt.xlabel(r'$x$')
plt.ylabel(r'$y$')


<class 'numpy.ndarray'>
Out[7]:
<matplotlib.text.Text at 0x7f8d3dc754a8>

Вопрос 2. Из скольки отрезков состоит график, изображающий предсказания построенного дерева на отрезке [-2, 2]?

  • 5
  • 6
  • 7
  • 8

2. Построение дерева решений для прогноза сердечно-сосудистых заболеваний

Считаем в DataFrame знакомый нам набор данных по сердечно-сосудистым заболеваниям.


In [3]:
df = pd.read_csv('../../mlcourse_open/data/mlbootcamp5_train.csv', 
                 index_col='id', sep=';')

In [4]:
df.head()


Out[4]:
age gender height weight ap_hi ap_lo cholesterol gluc smoke alco active cardio
id
0 18393 2 168 62.0 110 80 1 1 0 0 1 0
1 20228 1 156 85.0 140 90 3 1 0 0 1 1
2 18857 1 165 64.0 130 70 3 1 0 0 0 1
3 17623 2 169 82.0 150 100 1 1 0 0 1 1
4 17474 1 156 56.0 100 60 1 1 0 0 0 0

Сделайте небольшие преобразования признаков: постройте признак "возраст в годах", а также постройте по 3 бинарных признака на основе cholesterol и gluc, где они, соответственно, равны 1, 2 или 3. Эта техника называется dummy-кодированием или One Hot Encoding (OHE), удобней всего в данном случае использовать pandas.get_dummmies.


In [5]:
df['age_of_year'] = (df['age'] // 365.25).astype(int)

df.drop(['age'], axis=1, inplace=True)

df = pd.get_dummies(df, columns=['cholesterol', 'gluc'])

df.head()


Out[5]:
gender height weight ap_hi ap_lo smoke alco active cardio age_of_year cholesterol_1 cholesterol_2 cholesterol_3 gluc_1 gluc_2 gluc_3
id
0 2 168 62.0 110 80 0 0 1 0 50 1 0 0 1 0 0
1 1 156 85.0 140 90 0 0 1 1 55 0 0 1 1 0 0
2 1 165 64.0 130 70 0 0 0 1 51 0 0 1 1 0 0
3 2 169 82.0 150 100 0 0 1 1 48 1 0 0 1 0 0
4 1 156 56.0 100 60 0 0 0 0 47 1 0 0 1 0 0

Разбейте выборку на обучающую и отложенную (holdout) части в пропорции 7/3. Для этого используйте метод sklearn.model_selection.train_test_split, зафиксируйте у него random_state=17.


In [6]:
y = df['cardio'].astype(int)
X = df.drop('cardio', axis=1)

X.shape, y.shape


Out[6]:
((70000, 15), (70000,))

In [7]:
X_train, X_valid, y_train, y_valid = train_test_split(X, y, train_size=0.7, random_state=17)

In [8]:
X_train.shape, X_valid.shape, y_train.shape, y_valid.shape


Out[8]:
((49000, 15), (21000, 15), (49000,), (21000,))

Обучите на выборке (X_train, y_train) дерево решений с ограничением на максимальную глубину в 3. Зафиксируйте у дерева random_state=17. Визуализируйте дерево с помошью sklearn.tree.export_graphviz, dot и pydot. Пример дан в статье под спойлером "Код для отрисовки дерева". Обратите внимание, что команды в Jupyter notebook, начинающиеся с восклицательного знака – это терминальные команды (которые мы обычно запускаем в терминале/командной строке).


In [9]:
DT = DecisionTreeClassifier(max_depth=3, random_state=17)

In [10]:
DT.fit(X_train, y_train)


Out[10]:
DecisionTreeClassifier(class_weight=None, criterion='gini', max_depth=3,
            max_features=None, max_leaf_nodes=None,
            min_impurity_split=1e-07, min_samples_leaf=1,
            min_samples_split=2, min_weight_fraction_leaf=0.0,
            presort=False, random_state=17, splitter='best')

In [11]:
# используем .dot формат для визуализации дерева
export_graphviz(DT, feature_names=X.columns,
                out_file='../img/3.3.dot', filled=True)
!dot -Tpng ../img/3.3.dot -o ../img/3.3.png
!rm ../img/3.3.dot


/bin/sh: 1: dot: not found

Вопрос 3. Какие 3 признака задействуются при прогнозе в построенном дереве решений? (то есть эти три признака "можно найти в дереве")

  • weight, height, gluc=3
  • smoke, age, gluc=3
  • age, weight, chol=3
  • age, ap_hi, chol=3

Сделайте с помощью обученного дерева прогноз для отложенной выборки (X_valid, y_valid). Посчитайте долю верных ответов (accuracy).


In [12]:
DT.score(X_valid, y_valid)


Out[12]:
0.72128571428571431

Теперь на кросс-валидации по выборке (X_train, y_train) настройте глубину дерева, чтобы повысить качество модели. Используйте GridSearchCV, 5-кратную кросс-валидацию. Зафиксируйте у дерева random_state=17. Перебирайте параметр max_depth от 2 до 10.


In [13]:
tree_params = {'max_depth': list(range(2, 11))}

tree_grid = GridSearchCV(DT, tree_params, cv=5, n_jobs=-1)

In [14]:
%%time

tree_grid.fit(X_train, y_train)


CPU times: user 300 ms, sys: 70 ms, total: 370 ms
Wall time: 2.07 s
Out[14]:
GridSearchCV(cv=5, error_score='raise',
       estimator=DecisionTreeClassifier(class_weight=None, criterion='gini', max_depth=3,
            max_features=None, max_leaf_nodes=None,
            min_impurity_split=1e-07, min_samples_leaf=1,
            min_samples_split=2, min_weight_fraction_leaf=0.0,
            presort=False, random_state=17, splitter='best'),
       fit_params={}, iid=True, n_jobs=-1,
       param_grid={'max_depth': [2, 3, 4, 5, 6, 7, 8, 9, 10]},
       pre_dispatch='2*n_jobs', refit=True, return_train_score=True,
       scoring=None, verbose=0)

In [15]:
tree_grid.best_params_


Out[15]:
{'max_depth': 6}

Нарисуйте график того, как меняется средняя доля верных ответов на кросс-валидации в зависимости от значения max_depth.


In [16]:
max_deep = [x['max_depth'] for x in tree_grid.cv_results_['params']]
std_score = tree_grid.cv_results_['std_test_score']

plt.plot(max_deep, std_score)

plt.xlabel(r'max deep')
plt.ylabel(r'std score')


Out[16]:
<matplotlib.text.Text at 0x7f188c1c18d0>

Выведите лучшее значение max_depth, то есть такое, при котором среднее значение метрики качества на кросс-валидации максимально. Посчитайте также, какова теперь доля верных ответов на отложенной выборке. Все это можно сделать с помощью обученного экземпляра класса GridSearchCV.


In [18]:
acc1 = DT.score(X_valid, y_valid)
acc2 = tree_grid.best_estimator_.score(X_valid, y_valid)

(acc2 - acc1) / acc1 * 100


Out[18]:
0.60738099953786129

Вопрос 4. Имеется ли на кривой валидации по максимальной глубине дерева четкий пик, если перебирать max_depth от 2 до 10? Повысила ли настройка глубины дерева качество классификации (accuracy) более чем на 1% на отложенной выборке?

  • да, да
  • да, нет
  • нет, да
  • нет, нет

Обратимся опять (как и в 1 домашке) к картинке, демонстрирующей шкалу SCORE для расчёта риска смерти от сердечно-сосудистого заболевания в ближайшие 10 лет.

Создайте бинарные признаки, примерно соответствующие этой картинке:

  • $age \in [45,50), \ldots age \in [60,65) $ (4 признака)
  • верхнее артериальное давление: $ap\_hi \in [120,140), ap\_hi \in [140,160), ap\_hi \in [160,180),$ (3 признака)

Далее будем строить дерево решений с этим признаками, а также с признаками smoke, cholesterol и gender. Из признака cholesterol надо сделать 3 бинарных, соотв-х уникальным значениям признака ( cholesterol=1, cholesterol=2 и cholesterol=3), эта техника называется dummy-кодированием или One Hot Encoding (OHE). Признак gender надо перекодировать: значения 1 и 2 отобразить на 0 и 1. Признак лучше переименовать в male (0 – женщина, 1 – мужчина). В общем случае кодирование значений делает sklearn.preprocessing.LabelEncoder, но в данном случае легко обойтись и без него.

Итак, дерево решений строится на 12 бинарных признаках.

Постройте дерево решений с ограничением на максимальную глубину = 3 и обучите его на всей исходной обучающей выборке. Используйте DecisionTreeClassifier, на всякий случай зафикисровав random_state=17, остальные аргументы (помимо max_depth и random_state) оставьте по умолчанию.

Вопрос 5. Какой бинарный признак из 12 перечисленных оказался самым важным для обнаружения ССЗ, то есть поместился в вершину построенного дерева решений?

  • Верхнее артериальное давление от 160 до 180 (мм рт.ст.)
  • Пол мужской / женский
  • Верхнее артериальное давление от 140 до 160 (мм рт.ст.)
  • Возраст от 50 до 55 (лет)
  • Курит / не курит
  • Возраст от 60 до 65 (лет)

In [23]:
data = df.copy()

data["age_45-50"] = data["age_of_year"].apply(lambda x: 1 if x >= 45 and x < 50 else 0)
data["age_50-55"] = data["age_of_year"].apply(lambda x: 1 if x >= 50 and x < 55 else 0)
data["age_55-60"] = data["age_of_year"].apply(lambda x: 1 if x >= 55 and x < 60 else 0)
data["age_60-65"] = data["age_of_year"].apply(lambda x: 1 if x >= 60 and x < 65 else 0)

data["ap_hi_120-140"] = data["ap_hi"].apply(lambda x: 1 if x >= 120 and x < 140 else 0)
data["ap_hi_140-160"] = data["ap_hi"].apply(lambda x: 1 if x >= 140 and x < 160 else 0)
data["ap_hi_160-180"] = data["ap_hi"].apply(lambda x: 1 if x >= 160 and x < 180 else 0)

data["male"] = data["gender"].map({1: 0, 2: 1})
 
#df = pd.get_dummies(df, columns=['cholesterol'])

data.drop(["height", "weight", "ap_hi", "ap_lo", "alco", "active", 
    "ap_hi", "gender", "gluc_1", "gluc_2", "gluc_3", "age_of_year"], axis=1, inplace=True)

data.head()


Out[23]:
smoke cardio cholesterol_1 cholesterol_2 cholesterol_3 age_45-50 age_50-55 age_55-60 age_60-65 ap_hi_120-140 ap_hi_140-160 ap_hi_160-180 male
id
0 0 0 1 0 0 0 1 0 0 0 0 0 1
1 0 1 0 0 1 0 0 1 0 0 1 0 0
2 0 1 0 0 1 0 1 0 0 1 0 0 0
3 0 1 1 0 0 1 0 0 0 0 1 0 1
4 0 0 1 0 0 1 0 0 0 0 0 0 0

In [24]:
y = data['cardio'].astype(int)
X = data.drop('cardio', axis=1)

X.shape, y.shape


Out[24]:
((70000, 12), (70000,))

In [25]:
X_train, X_valid, y_train, y_valid = train_test_split(X, y, train_size=0.7, random_state=17)

In [26]:
X_train.shape, X_valid.shape, y_train.shape, y_valid.shape


Out[26]:
((49000, 12), (21000, 12), (49000,), (21000,))

In [27]:
tree = DecisionTreeClassifier(max_depth=3, random_state=17)

In [28]:
tree.fit(X_train, y_train)


Out[28]:
DecisionTreeClassifier(class_weight=None, criterion='gini', max_depth=3,
            max_features=None, max_leaf_nodes=None,
            min_impurity_split=1e-07, min_samples_leaf=1,
            min_samples_split=2, min_weight_fraction_leaf=0.0,
            presort=False, random_state=17, splitter='best')

In [29]:
# используем .dot формат для визуализации дерева
export_graphviz(tree, feature_names=X.columns,
                out_file='../img/3.5.dot', filled=True)
!dot -Tpng ../img/3.5.dot -o ../img/3.5.png
!rm ../img/3.5.dot